library(sf)
library(dplyr)
library(ggplot2)
library(scales)
library(ggmap)
library(leaflet)

Prepare regions shapefile

ak_regions <- read_sf("data/shapefiles/ak_regions_simp.shp")

st_crs(ak_regions)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

class(ak_regions)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

#plot(ak_regions)

Use different projection with Alaska Albers

ak_regions_3338 <- ak_regions %>%
  st_transform(crs = 3338)

st_crs(ak_regions_3338)

## Coordinate Reference System:
##   EPSG: 3338 
##   proj4string: "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

plot(ak_regions_3338)

summary(ak_regions_3338)

##    region_id     region            mgmt_area          geometry 
##  Min.   : 1   Length:13          Min.   :1   MULTIPOLYGON :13  
##  1st Qu.: 4   Class :character   1st Qu.:2   epsg:3338    : 0  
##  Median : 7   Mode  :character   Median :3   +proj=aea ...: 0  
##  Mean   : 7                      Mean   :3                     
##  3rd Qu.:10                      3rd Qu.:4                     
##  Max.   :13                      Max.   :4

Prepare population data

pop <- read.csv("data/shapefiles/alaska_population.csv", stringsAsFactors = FALSE)
head(pop)

##   year     city      lat       lng population
## 1 2015     Adak 51.88000 -176.6581        122
## 2 2015   Akhiok 56.94556 -154.1703         84
## 3 2015 Akiachak 60.90944 -161.4314        562
## 4 2015    Akiak 60.91222 -161.2139        399
## 5 2015   Akutan 54.13556 -165.7731        899
## 6 2015 Alakanuk 62.68889 -164.6153        777

Coerce the sf object using st_as_sf

pop_4326 <- st_as_sf(pop, coords = c("lng", "lat"), crs = 4326, remove = F)
st_crs(pop_4326)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

head(pop_4326)

## Simple feature collection with 6 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -176.6581 ymin: 51.88 xmax: -154.1703 ymax: 62.68889
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   year     city      lat       lng population                   geometry
## 1 2015     Adak 51.88000 -176.6581        122    POINT (-176.6581 51.88)
## 2 2015   Akhiok 56.94556 -154.1703         84 POINT (-154.1703 56.94556)
## 3 2015 Akiachak 60.90944 -161.4314        562 POINT (-161.4314 60.90944)
## 4 2015    Akiak 60.91222 -161.2139        399 POINT (-161.2139 60.91222)
## 5 2015   Akutan 54.13556 -165.7731        899 POINT (-165.7731 54.13556)
## 6 2015 Alakanuk 62.68889 -164.6153        777 POINT (-164.6153 62.68889)

Change coord reference system

pop_3338 <- pop_4326 %>%
  st_transform(crs = 3338)

Do the spatial join

pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)
head(pop_joined)

## Simple feature collection with 6 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -1537925 ymin: 472627.8 xmax: -10340.71 ymax: 1456223
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##   year     city      lat       lng population region_id           region
## 1 2015     Adak 51.88000 -176.6581        122         1 Aleutian Islands
## 2 2015   Akhiok 56.94556 -154.1703         84         6           Kodiak
## 3 2015 Akiachak 60.90944 -161.4314        562         8        Kuskokwim
## 4 2015    Akiak 60.91222 -161.2139        399         8        Kuskokwim
## 5 2015   Akutan 54.13556 -165.7731        899         1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153        777        13            Yukon
##   mgmt_area                   geometry
## 1         3  POINT (-1537925 472627.8)
## 2         3 POINT (-10340.71 770998.4)
## 3         4  POINT (-400885.5 1236460)
## 4         4  POINT (-389165.7 1235475)
## 5         3 POINT (-766425.7 526057.8)
## 6         4  POINT (-539724.9 1456223)

Calculate population by region

pop_region <- pop_joined %>%
  as.data.frame() %>%
  group_by(region) %>%
  summarise(total_pop = sum(population))

head(pop_region)

## # A tibble: 6 x 2
##   region           total_pop
##   <chr>                <int>
## 1 Aleutian Islands      8840
## 2 Arctic                8419
## 3 Bristol Bay           6947
## 4 Chignik                311
## 5 Cook Inlet          408254
## 6 Copper River          2294

pop_region_3338 <- left_join(ak_regions_3338, pop_region)

## Joining, by = "region"

head(pop_region_3338)

## Simple feature collection with 6 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -2175328 ymin: 405653.9 xmax: 773820 ymax: 2383770
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 6 x 5
##   region_id region   mgmt_area                           geometry total_pop
##       <int> <chr>        <dbl>                 <MULTIPOLYGON [m]>     <int>
## 1         1 Aleutia…         3 (((-1156666 420855.1, -1159837 41…      8840
## 2         2 Arctic           4 (((571289.9 2143072, 569941.5 214…      8419
## 3         3 Bristol…         3 (((-339688.6 973904.9, -339302 97…      6947
## 4         4 Chignik          3 (((-114381.9 649966.8, -112866.8 …       311
## 5         5 Copper …         2 (((561012 1148301, 559393.7 11481…      2294
## 6         6 Kodiak           3 (((115112.5 983293, 113051.3 9828…      8126

#plot(pop_region_3338)

pop_mgmt <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarise(total_pop = sum(total_pop), do_union = F)

plot(pop_mgmt["total_pop"])

Make maps

rivers_3338 <- read_sf("data/shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

head(rivers_3338)

## Simple feature collection with 6 features and 3 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: -622169.8 ymin: 557375.7 xmax: 738652.5 ymax: 2316718
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## # A tibble: 6 x 4
##   StrOrder region          n                                       geometry
##      <int> <chr>       <int>                          <MULTILINESTRING [m]>
## 1        3 Aleutian I…   316 ((-616509.9 557976.9, -616240.8 557375.7), (-…
## 2        3 Arctic      10869 ((-20959.19 2009003, -20804.35 2009054), (-20…
## 3        3 Bristol Bay  4884 ((-231927.5 794449.9, -232004.6 794730.2), (-…
## 4        3 Chignik       141 ((-302874.8 715936.3, -302829.4 715920.8), (-…
## 5        3 Cook Inlet   3447 ((-11441.06 987666.2, -11772.61 988142.6), (-…
## 6        3 Copper Riv…  1784 ((540428.6 1180052, 539410.1 1176802), (54043…

ggplot() +
  geom_sf(data = pop_region_3338, mapping = aes(fill = total_pop)) +
  geom_sf(data = pop_3338, mapping = aes(), size = .5) +
  geom_sf(data = rivers_3338, mapping = aes(size = StrOrder), color = "Black") +
  scale_size(range = c(.01, .2), guide = F) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma)

Incorporate base maps using ggmap

pop_3857 <- pop_3338 %>%
  st_transform(crs = 3857)

# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}

bbox <- c(-170, 52, -130, 64)
ak_map <- get_stamenmap(bbox, zoom = 4)

## Source : http://tile.stamen.com/terrain/4/0/4.png

## Source : http://tile.stamen.com/terrain/4/1/4.png

## Source : http://tile.stamen.com/terrain/4/2/4.png

## Source : http://tile.stamen.com/terrain/4/0/5.png

## Source : http://tile.stamen.com/terrain/4/1/5.png

## Source : http://tile.stamen.com/terrain/4/2/5.png

ak_map_3857 <- ggmap_bbox_to_3857(ak_map)

ggmap(ak_map_3857) +
  geom_sf(data = pop_3857, aes(color = population), inherit.aes = FALSE) +
  scale_color_continuous(low = "khaki", high = "firebrick", labels = comma)

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Leaflet

epsg3338 <- leaflet::leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:3338",
  proj4def =  "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(16:7))

pop_region_4326 <- pop_region_3338 %>% 
  st_transform(4326)

leaflet(options = leafletOptions(crs = epsg3338)) %>%
  addPolygons(data = pop_region_4326, fill = "gray", weight = 1)

pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")

m

pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1) %>% 
        addCircleMarkers(data = pop_4326,
                         lat = ~lat,
                         lng = ~lng,
                         radius = ~log(population/500), # arbitrary scaling
                         fillColor = "gray",
                         fillOpacity = 1,
                         weight = 0.25,
                         color = "black",
                         label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")

m

Spatial Analysis

Ben Gerber

11/8/2019

Prepare regions shapefile

Prepare population data

Calculate population by region

Make maps

Incorporate base maps using ggmap

Leaflet